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1.  Introduction 


Viscous  shock-on-shock  interactions  occur  on  typical  aerospace  configurations  such  as  the  proposed 
NASP  vehicle  and  the  space  shuttle.  The  forebodies  of  proposed  hypersonic  aircraft  (Fig.  1.1) 
form  ramp-like  structures  designed  to  compress  the  incoming  air  with  oblique  shocks  and  there!'  ,re 
act  as  the  compressor  system  for  the  inlet.  For  optimum  mass  flow  through  the  inlet,  tlie.se 
compression  system  shocks,  which  may  form  a  relatively  strong  oblique  shock  in  conjunction  with 
the  vehicle  bow  shock,  should  be  positioned  to  converge  on  the  inlet  cowl  leading  edge  [1]  wlnue 
they  interact  with  the  bow  shock  produced  by  the  cowl  lip.  Viscous  hypersonic  shock-on-shock 
interactions  (often  denoted  “interfering”  flows)  can  significantly  affect  the  performance  of  the  inlet 
through  the  creation  of  anomalous  pressure  and  heat  transfer  peaks  on  the  cowl  leading  edge. 

The  ability  to  accurately  and  efficiently  predict  flow  characteristics  of  interfering  flows  is 
critical  for  successful  design  [2].  This  effort  has  therefore  received  significant  impetus  in  recent 
years  under  the  auspices  of  the  NASP  Project.  A  recent  report,  Hypersonic  Technology  for 


ta*y  nppii^ctti oiio  ^j,  emphasizes  the  importance  oc  tuc  study  oi  ouch  interactions; 


The  most  intense  local  heating  rates  on  vehicles  such  as  the  projected  NASP  research 
vehicle  are  expected  to  be  on  cowl  lips,  caused  by  shock-on-shock  heating. 


Edney  [4]  classified  shock  interference  patterns  into  six  types.  The  actual  pattern  established 
in  a  particular  case  depends  primarily  upon  the  strength  and  location  of  the  impinging  shock, 
the  characteristics  of  the  incoming  flow  and  the  cowl  shape.  Efforts  at  understanding  interfering 
flows  have  focused  on  experimental  [4,  5,  6,  7],  semiempirical  [8],  theoretical  [9,  10]  and,  in  recent 
years,  numerical  simulations  [1,  2,  11,  12,  13]. 

Most  studies  so  far  have  focused  on  Type  III  and  Type  IV  interactions  which  result  in  ihe 
highest  mechanical  and  thermal  loading  on  the  cowl  lip  (Fig.  1.2).  These  types  of  interactions 
occur  when  the  impinging  shock  intersects  the  subsonic  portion  of  the  cowl  bow  shock  —  the  cowl 
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is  typically  replaced  by  a  cylinder  to  simplify  analysis.  A  supersonic  viscous  shear  layer  emanates 
from  the  point  of  impingement  and,  depending  upon  the  angle  between  the  shear  layer  and  the 
surface  tangent  at  the  point  where  the  shear  layer  hits  the  surface,  this  layer  either  merges  with 
the  boundary  layers  where  it  strikes  the  body  and  becomes  an  attached  subsonic  layer  (Type  111 ) 
or.  the  shear  layer  forms  a  “jet”  which  strikes  the  body  after  the  formation  of  a  terminating  jot 
bow  shock  (Type  IV)  [1,  10].  The  process  of  jet  impingement  (Type  /V  )  results  in  maximum 
augmentation  of  peak  heat-transfer  and  surface  pressure. 

Initial  theoretical/numerical  efforts,  relied  upon  semiempirical  approaches  requiring  specifica¬ 
tion  of  several  parameters  such  as  shock  standoff  distance.  Hains  and  Keyes  [14],  for  example, 
utilized  several  assumptions  to  develop  empirical  correlations  for  peak  heat  transfer  amplifica¬ 
tion.  Based  on  a  study  of  flow  past  simple  body  shapes  at  different  Mach  numbers  (0.0  and  20.2) 
and  specific  heat  ratios  (1.2  to  1.6),  they  concluded  that  effect  of  an  impinging  shock  is  most 
drastic  for  Type  IV  interactions  with  amplification  up  to  17  times  (over  the  noninterfering  case) 
•'or  heat-transfer  rate  and  8  times  for  peak  pressure.  Other  notable  efforts  include  the  work  of 
Crawford  [9]  who  developed  a  graphical  method  of  pattern  prediction  and  Biamiette  [10]  who 
introduced  further  simplifications  to  Crawford’s  method. 

The  advent  of  more  sophisticated  methods  and  high-speed  computers  saw  significant  advances 
in  the  use  of  computational  methods.  Tannehill  et  al.  [11]  solved  the  full  2-D  Navier-Stokes  equa¬ 
tions  with  shock  capturing  as  well  as  shock  fitting  methods  for  low  and  high  Reynolds  Numbers 
respectively  for  Type  III  and  IV  flows  around  a  cylinder  at  Mach  4.6.  Although  their  explicit 
algorithm  suffered  significant  step  size  limitations,  they  successfully  validated  their  approach  by 
comparison  with  the  results  of  Beckwith  and  Cohen  [15]  and  Edney  [4].  White  and  Rhie  [1]  com¬ 
puted  blunt  cowl  flows  with  (Mach  6)  and  without  (Mach  15)  shock  impingement  with  a  pressure 
based  implicit  finite  volume  scheme.  Their  results  compared  well  with  the  experiments  of  Craig 
and  Ortwerth  [6]  and  Tannehill  et  al  [11]. 

The  complexity  of  the  wave  structure  of  shock-on-shock  interactions  of  the  type  investigated 
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in  this  research  all  but  obviates  the  use  of  shock-fitting  methods.  Modem  analyses  utilize  s ! .  ■  <  k 
capturing  methods  almost  exclusively.  The  first  few  popular  shock-capturing  schemes  such  ns 
the  centered  MacCormack  scheme  [16]  utilized  by  Tannehill  et  al  or  the  Jameson  scheme  '17] 
require  the  addition  of  explicit  damping  terms  to  avoid  overshoots  and  oscillations  leading  otom 
to  extremely  smeared  shocks  for  hypersonic  flows. 

The  underlying  philosophy  of  more  recent  shock-capturing  methods  is  either  mathemetical 
(for  example  TVD  schemes  [18])  or  physical  (for  example  flux-splitting  [19,  20]).  Many  mod¬ 
ern  schemes  employ  upwinding  to  obtain  algorithms  possessing  better  dissipation  characteristics, 
higher  stability  bounds  and  increased  numerical  efficiency  [21,  22].  A  brief  description  of  modern 
schemes  utilized  in  the  recent  past  for  the  types  of  flows  under  investigation  is  presented  below. 

TVD  schemes,  upwind  biased  and  symmetric,  have  gained  rapid  popularity  in  recent  years  [23] 
and  have  been  applied  for  several  flows  including  shock  wave  diffraction,  flows  past  airfoils,  com¬ 
plex  vehicle  shapes  [24],  boundary  layers  and  shock  boundary  layer  interactions  [25].  Klopfer  .  id 
Yee  [2]  computed  noninterfering  and  interfering  patterns  of  all  types  identified  by  Edney  [  I  for 
Mach  6.  8  and  15  flows.  Utilizing  the  thin-layer  Navier  Stokes  equations,  they  suggested  grid  res¬ 
olution  criteria  for  noninterfering  flows  and  concluded  that  their  scheme  performed  relatively  well 
for  2-D  flows  though  modest  discrepancies  existed  in  comparison  with  experimental  heat  transfer 
data. 

The  popular  van  Leer  [20]  and  Steger-Warming  [19]  flux  split  schemes  have  been  applied 
to  several  inviscid  and  viscous  flows  in  explicit  and  implicit  formulations.  Applications  with  van 
Leer’s  scheme  include  inviscid  subsonic  and  transonic  flows  over  airfoils  [26],  viscous  shock-induced 
separated  flows  [27]  and  flows  over  delta  wings  [28].  Recently,  Moon  and  Holt  [12]  reported  a  com¬ 
putation  with  van  Leer  splitting  for  inviscid  fluxes  and  a  centered  scheme  for  viscous  terms  for  a 
Type  III+  interaction  (Mach  8)  with  and  without  a  turbulence  model.  Their  results  indicate  sig¬ 
nificant  overprediction  of  surface  peak  pressures  for  the  laminar  case  and  under  prediction  of  peak 
heat  transfer  for  both  cases.  In  addition,  the  turbulent  computation  displavs  spatial  oscillations 
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in  the  region  of  maximum  heat  transfer. 

Reported  applications  of  the  Steger- Warming  algorithm  include  inviscid  flows  past  airfoils  [19] 
and  cylinders  [29],  The  numerical  algorithm  employed  in  the  current  effort  is  a  modified  Steger- 
U arming  scheme  proposed  by  MacCormack  and  Candler  [30]  (Section  3).  This  scheme  has  previ¬ 
ously  been  applied  for  flat  plate  boundary  layer  flows,  flows  past  compression  ramps,  blunt  body 
(lows  [30]  and  for  viscous  real  gas  flows  past  sphere-cones  [31].  A  comparison  of  the  van  Leer  and 
the  original  Steger- Warming  scheme  for  some  Euler  flows  may  be  found  in  Anderson  et  al.  [32]. 

Recent  research  in  the  computation  of  hypersonic  flows  has  also  focused  on  unstructured  grid 
methods  which  are  typically  well-suited  to  adaptive  techniques.  Thareja  et  al.  [13]  presented  a 
2-1)  upwind  finite  element  technique  using  cell  centered  quantities  and  implicit  and/or  explicit 
time  marching  with  adaptive  unstructured  triangular  grids.  They  implemented  a  first  order  basic 
and  a  higher  order  flux  corrected  scheme  [33]  with  an  essentially  point  Gauss-Seidel  implicit 
algorithm.  The  shocks  obtained  with  this  method  are  crisp  though  some  smear  is  evident  in  the 
-hear  layers.  Surface  quantities  including  pressure  and  heat  transfer  rates  compare  very  well  with 
mi  her  computational  schemes  and  experimental  data. 

This  report  first  outlines  the  objectives  of  this  research  effort.  This  is  followed  by  a  relatively 
detailed  exposition  of  the  numerical  algorithm  in  so  far  as  the  inviscid  flux  evaluation  is  con¬ 
cerned.  For  the  purposes  of  brevity,  the  Gauss-Seidel  line  relaxation  algorithm  employed  for  time 
advancement  is  only  summarized.  The  computations  described  are  then  classified  prior  to  the 
discussion  of  the  results. 
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TYPICAL  CONFIGURATION 


TYPICAL  OVERALL  FLOW  FIELD 


Figure  1.1:  Schematic  of  NAS P  vehicle.  From  [2] 


BOW  SHOCK 


Figure  1.2:  Type  III  and  IV  shock  interference  patterns.  From  [7] 


2.  Objective 


The  computation  of  shock-on-shock  interactions  is  a  challenging  task  requiring  the  resolution  of 
intense  shock  waves  and  slip  streams  and  their  interactions  with  each  other  and  the  boundary 
layer.  As  mentioned  earlier,  interfering  viscous  hypersonic  flows  past  cylinders  have  been  com¬ 
puted  previously  with  TVD  schemes  [2],  van  Leer  splitting  [12]  and  finite  element  flux  corrected 
transport  methods  [13].  Each  of  these  methods  possesses  some  degree  of  built  in  or  explicitly 
added  damping  to  stabilize  the  calculation  in  regions  of  high  gradients  such  as  shock  waves.  A 
brief  description  of  the  damping  properties  of  several  popular  schemes  is  presented  in  Ref.  [34]. 
Several  modern  schemes  were  in  fact  developed  primarily  to  handle  inviscid  dominated  phenomena 
and  focused  almost  exclusively  on  wave  propagation  by  considering  only  the  inviscid  terms  in  the 
governing  equations.  The  introduction  of  viscous  terms  (typically  through  central  differencing) 
introduces  contradictory  requirements  on  the  numerical  scheme  i.e.,  shock-capturing  requires  a  fi¬ 
nite  amount  of  numerical  dissipation  which  must  not  however,  overwhelm  the  physical  dissipation 
in  the  boundary  layer. 

The  fundamental  objective  of  this  study  is  to  utilize  the  modified  Steger- Warming  flux  split 
scheme  of  MacCormack  and  Candler  [30]  to  examine  viscous  shock-on-shock  impingement  flow- 
fields.  A  salient  objective  is  to  examine  the  characteristics  of  the  scheme  with  grid  resolution 
studies  and  to  compare  this  scheme  with  others.  The  choice  of  scheme  is  motivated  by  the  fol¬ 
lowing  factors.  1)  The  modifications  specifically  address  the  issue  of  excessive  dissipation  in  the 
boundary  layer.  Such  undesirable  dissipation  can  significantly  degrade  the  evaluation  of  surface 
quantities.  2)  The  algorithm  has  been  extended  to  include  non-equilibrium  real  gas  and  success¬ 
fully  applied  for  hypersonic  blunt  body  flows  [35].  It  is  an  attractive  candidate  for  further  work  on 
complex  hypersonic  3-D  configurations  and  thus,  its  capabilities  and  limitations  warrant  further 
investigation  especially  under  complex  flow  conditions. 

The  configurations  examined  in  this  paper  are  necessarily  dictated  by  the  availability  of  ex- 
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isting  experimental  and  numerical  investigations.  For  the  purposes  of  comparison  with  other 
schemes,  the  computations  of  Moon  and  Holt  [12]  (Type  ///+)  and  Thareja  et  al.  [13]  (Type  IV) 
are  utilized  as  described  later. 
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3.  Theoretical  Model 


The  2-D  Navier  Stokes  equations  in  strong  conservation  form  are  solved  in  the  transformed  (£.  ?/) 


coordinates: 


8U_  d£  dG_n 

dt  d£  drj 


where  U,  F  and  G  are  vectors  as  defined  below  and  £(x,y)  and  77(1,3/)  are  the  transformed 
variables  (with  d£/dx  denoted  by  etc.).  The  solution  vector  U  is  written  as: 


*-7 


and  the  flux  vectors  F  and  G  are: 


1  I  puU  +  £x(P  Txx)  £yTXy 
F  =  T 

pvU  +  £y(p  ~  Tyy)  ~  txTxy 

.  (j Oe  +  p)U  +  £Xf3X  -  ZyPy 


1  puV  +  T] x(p  Txx )  Vy^xy 

G=l\  I ’ 

PVV  +  T]y(p  ~  Tyy)  -  TJ^y 

k  (PC  +  P)1 * * * V  +  VxPx  ~  VyPy  , 

and,  in  Eqns.  3.3  and  3.4,  the  contravariant  velocity  vectors  U  and  V  are: 


U  —  £XU  +  £yV 

V  =  t]xU  +  T]yV 


and  J  is  the  Jacobian  of  the  transformation,  with  subscripts  indicating  derivatives  (£x  =  dti/dx): 

J~l  =  Xtyv  -  Xr,yz  (3.7) 


The  components  of  the  Cartesian  stress  tensor  r  may  be  written  as, 

,  _  du 

Txx  =  Ay.i)  +  2/i— 
ox 

.du  dv 
r*»  =  + 


‘yy 


,  „  n  dv 

=  +  2(i— 

dy 


where  A  =  — 2/3/x,  and 


Px  —  Qx  UTXX  VTXy 
Py  ~  Qy  ~  UTxy  ~  VTyy 
The  components  of  the  heat  flux  vector  are 

H  dt{ 


Qx  = 
Qy  = 


^  Pr  dx 
H  dei 
Pr  dy 


(3.8) 

(3.9) 

(3.10) 

(3.11) 

(3.12) 

(3.13) 

(3.14) 


The  Cartesian  velocity  components  in  the  ( x,y )  coordinates  are  denoted  by  (u,  v).  For  the 
configurations  under  consideration,  the  7/  coordinate  lines  are  in  the  generally  radial  direction 
while  the  £  lines  are  circumferential.  The  density  p,  static  pressure  p  and  static  temperature  T 
are  related  through  the  equation  of  state  p  =  pRT  where  R  is  the  gas  constant.  The  total  energy 
per  unit  mass  is  given  by  e  =  +  0.5 (u2  +  v2  -f  w2),  where  the  internal  energy  per  unit  mass  e, 

is  equal  to  CVT.  The  molecular  dynamic  viscosity  p  is  given  by  the  Sutherland’s  law: 


P  =  Ci- 


T2 


r  +  c2  (3-15» 

where  C\  and  C2  are  constants  (=  1.45  x  10 ~6kg/(m  s-/K)  and  110. 4A'  respectively).  The 
molecular  Prandtl  number  Pr  =  Cpp / k  is  taken  to  be  0.73  (air). 
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The  above  equations  are  solved  in  discretized  form  with  the  cell-centered  finite  volume  Gauss- 
Seidel  line  relaxation  numerical  algorithm  described  in  [36]  utilizing  the  approach  described  in  [30] 
to  evaluate  the  inviscid  fluxes.  A  brief  development  of  the  algorithm  is  presented  for  completeness 
with  reference  to  the  evaluation  of  the  flux  at  the  cell  surface  j  +  \.  The  flux  computation  proceeds 
separately  for  the  viscous  and  inviscid  parts.  At  any  time  level,  suppressing  the  overbars, 

=  G/J+S  +  Gv.,+}  (3-16) 

where  the  subscripts  I  and  V  denote  the  inviscid  and  viscous  components  respectively.  Since  the 
inviscid  fluxes  are  homogeneous,  following  the  procedure  of  Steger  and  Warming  [19],  they  may 
be  split  into  subvectors  possessing  advantageous  eigenvalue  properties: 

G?  =  BnlP  =  {Bl  +  Bl)Tr  (3.17) 

Utilizing  the  hyperbolic  nature  of  the  inviscid  fluxes,  the  Jacobian  B  is  diagonalized  with  a 
similarity  transformation: 

B  =  Q-1\Q  =  Q-'(A+  +  A~)Q  =  B+  +  B _  (3.18) 

where  A  is  a  diagonal  matrix  consisting  of  the  eigenvalues  of  B  and  A+  and  A-  denote  the  splitting 
of  the  eigenvalues  into  positive  and  negative  components.  For  simplicity  of  evaluation,  Q  may 
further  be  written  as: 

Q  =  CRS  (3.19) 

where  R  is  a  rotation  matrix,  S  =  V  is  the  vector  of  non-conserved  variables  {p,  u,v,p)  and 
C  diagonalizes  the  flux  vector  G  written  in  terms  of  V .  At  a  face  j  +  1/2  therefore,  the  flux 
according  to  the  method  of  Steger  and  Warming  may  be  written  as: 

GIj+h  =  B+)Ut<J  +  (3.20) 

The  quantities  U  are  obtained  at  the  cell  interfaces  by  extrapolating  the  conserved  variables  to  the 
cell  surface  with  the  MUSCL  approach  of  Anderson  et  al.  [32]  as  described  below.  The  formula 
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of  Eqn.  3.20  exhibits  two  major  difficulties.  First,  it  may  introduce  discontinuities  in  the  solution 
at  sonic  and  stagnation  points  where  the  eigenvalues  change  sign  and  some  form  of  eigenvalue 
smoothing  may  be  necessary.  For  one  of  the  computations  reported  in  this  report,  smoothing  pro¬ 
portional  to  the  eigenvalue  gradients  (cf.  [32])  is  employed  (Section  6.2).  A  more  serious  problem 
introduced  by  the  evaluation  of  Eqn.  3.20  is  the  problem  of  excessive  numerical  damping  intro¬ 
duced  in  the  boundary  layers.  This  damping  significantly  deterioriates  the  evaluation  of  surface 
quantities  of  engineering  interest.  MacCormack  and  Candler  [30]  considered  the  determination  of 
the  normal  flux  near  a  surface  under  boundary  layer  conditions.  With  direct  algebraic  manipula¬ 
tion,  they  proved  that  the  Steger  Warming  procedure  introduces  artifical  tangential  momentum 
exchange  between  adjacent  cells  in  the  boundary  layer  solely  due  to  the  splitting  of  the  inviscid 
fluxes.  This  diffusive  term  is  proportional  to  the  quantity: 

2^Pi,j  ~  u*,j)  (3.21) 

Since  it  is  of  order  Ay,  it  can  obtain  unacceptably  high  values  in  the  boundary  layer.  In  addition, 
there  is  also  a  large  numerical  exchange  of  kinetic  energy  of  the  order: 

~  ai,j)  (3.22) 

between  adjacent  points  in  the  boundary  layer  where  a  is  the  kinetic  energy.  They  recommended 
that  both  J9+ ai,d  be  evaluated  at  the  same  point  (either  j  or  j  +  1).  In  the  present 

calculations,  following  this  recommendation: 


n  odd 
n  even 


(3.23) 


where  n  is  the  current  iteration  number.  Further  investigation  showed  however,  that  this  mod¬ 
ification  introduced  a  second  order  fictitious  pressure  gradient  across  the  boundary  layer.  This 
error  comprised  the  terms: 
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and  may  therefore  be  expected  to  dominate  only  in  the  close  vicinity  of  the  boundary  where 


high  velocity  gradients  exist.  MacCormack  and  Candler  corrected  this  gradient  at  the  expense  of 
reintroducing  some  of  the  diffusiveness  of  the  Steger  Warming  scheme  by  further  modifying  the 
component  S  for  the  matrix  Q  in  Eqn.  3.19  so  that  its  last  row  is  evaluated  as  in  the  original 
formulation.  In  this  work,  this  last  row  ( LR )  correction  is  applied  in  a  linear  fashion  starting  at 
the  nominal  edge  of  the  boundary  layer  obtaining  the  full  correction  at  the  surface. 

The  modifications  described  above  were  developed  for  application  only  in  the  boundary  layers 
and  in  fact  lead  to  instability  when  applied  in  regions  of  discontinuities  such  as  shock  waves.  As 
a  result,  it  is  necessary  to  revert  to  the  original  Steger  Warming  scheme  near  shocks.  Further,  in 
order  to  retain  the  monotonic  nature  of  the  solution,  it  is  also  necessary  to  revert  to  first  order 
accuracy  at  such  regions.  The  order  of  the  solution  is  determined  by  appropriate  extrapolation 
of  the  conserved  variables  as  in  the  MUSCL  approach  which  is  known  to  be  superior  to  the  flux- 
differencing  approach  in  which  the  split  fluxes  are  first  evaluated  at  the  nodal  points  and  t  hen 
extrapolated  to  the  cell  surfaces  [32], 


V]+i  =  l>i+*] 


(Vj  -  V 


=  -*u- 


+  (c,+2  -  0i+,) 


(3.25) 


(3.26) 


where  the  superscripts  -  and  +  indicate  states  to  the  left  and  right  respectively  of  the  interface 


j  +  2- 


Formal  first  and  second  order  accuracy  is  obtained  by  choosing  the  term  equal  to  zero  and 
one  respectively.  A  pressure  switch  given  by  (for  the  j  +  ^  face): 


Pi -Pi -i 

min{Pj,Pj- 1) 


otherwise 


(3.27) 


is  utilized  to  determine  regions  where  the  scheme  must  revert  to  the  first  order  accurate  Steger 
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Warming  method  for  monotonicity  and  stability.  The  value  of  A  is  typically  fixed  at  0.25.  Other 
strategies  are  described  in  [23,  32]. 

The  governing  equations  (Eqn.  3.1)  are  integrated  with  a  residual  driven  line  Gauss-Seidel 
relaxation  scheme  as  described  in  the  works  of  MacCormack  [36]  and  Candler  [35].  With  first- 
order  backward  Euler  time  discretization  and  linearization  of  the  fluxes  in  time,  the  discretized 
equation  may  be  written  in  the  form: 

{NUMERICS}6U  =  PHYSICS  (3.28) 

where  PHYSICS  represents  the  residual  and  NUMERICS  contains  the  driving  terms  and  <5(7 
represents  the  change  in  the  solution  vector  at  each  time  step.  The  full  Navier-Stokes  equations  are 
utilized  in  computing  the  residual.  The  viscous  terms  in  the  residual  are  evaluated  with  the  second 
order  (in  space  and  time)  predictor-corrector  algorithm  of  MacCormack[16].  One  advantage  of 
the  above  methodology  is  that  since  the  NU M ERICS  portion  of  the  code  is  simply  a  driver  for 
the  solution,  approximations  may  be  utilized.  In  fact,  Liou  and  van  Leer  [37]  point  out  that  even 
if  the  PHYSICS  term  is  evaluated  with  other  methods  (such  as  van  Leer’s  splitting  or  Roe’s 
upwinding),  the  use  of  Steger  Warming  splitting  in  evaluating  the  NU M ERICS  portion  leads 
to  robust  codes  for  the  Newton  linearization  procedure.  In  this  research,  the  viscous  Jacobians  in 
the  driving  NUM ERICS  portion  are  computed  with  the  thin  layer  approximation. 

When  solved  with  the  Gauss-Seidel  line  relaxation  procedure,  Eqn.  3.28  represents  a  block 
tridiagonal  system  displaying  strong  diagonal  dominance  when  the  above  flux-splitting  is  em¬ 
ployed.  For  transonic  and  supersonic  flows,  line  relaxation  methods  are  superior  to  approximate 
factorization  methods  in  convergence  rate  which  amply  compensates  for  the  larger  computation 
required  per  iteration  [38].  The  line  Gauss-Seidel  algorithm  is  also  unconditionally  stable  in 
the  linear  analysis  and  is  known  to  be  relatively  insensitive  to  the  choice  of  time  increment  per 
iteration. 
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4.  Classification  of  Computations 


The  geometrical  complexity  of  a  3-D  cowl  is  separated  from  the  physics  of  the  flow  features  by 
considering  a  simplified  geometry  as  shown  schematically  in  Fig  4.1.  The  cowl  is  replaced  by  a 
cylinder  in  an  oncoming  flow  containing  the  impinging  shock.  The  length  dimensions  and  flow 
parameters  chosen  (described  below)  are  determined  by  availability  of  experimental  data  and 
existing  numerical  studies  for  comparison  purposes.  The  characteristics  of  the  impinging  shock 
are  uniquely  specified  by  the  flow  deflection  angle  and  a  reference  point  in  the  oncoming  flow. 
Two  distinct  configurations  are  computed  with  the  main  distinction  being  in  the  shock  strength 
and  impingement  location. 

The  first  set  of  computations  corresponds  to  the  experimental  configuration  of  Wieting  and 
Holden  [5]  and  computed  with  an  adaptive  upwind  finite  element  technique  by  Thareja  et  al.  [13]. 
This  configuration,  representing  the  computation  described  in  most  detail  by  Thareja  et  al..  will 
be  denoted  henceforth  as  Config.  A.  The  second  set  of  computations,  denoted  Config.  B  is  based 
on  the  experiments  of  Wieting  [39]  and  computed  with  the  van  Leer  flux-split  scheme  by  Moon 
and  Holt  [12].  The  flow  parameters  and  geometries  are  described  in  Table  4.1.  Note:  For  ease  of 
comparison  with  previous  work,  results  for  configuration  A  are  shown  inverted  about  the  cylinder 
horizontal  center-line. 

For  each  configuration,  three  grids  providing  varying  degrees  of  resolution  in  both  the  f  and 
rj  directions  are  utilized.  The  inflow  boundary  [AB  in  Fig.  4.1]  is  assumed  to  be  a  sequence  of 
piecewise  smooth  spirals  bounding  the  domain  of  interest.  For  Config.  A,  the  three  grids  represent 
increased  concentration  throughout  the  domain  while  for  Config.  B,  the  finest  grid  concentrates 
more  points  only  in  the  region  of  interest.  The  cases  are  denoted  Al,  A2  and  A3  representing 
the  coarsest,  medium  and  finest  grids  respectively  with  similar  nomenclature  for  Config.  B.  In 
the  radial  direction  the  grids  are  stretched  algebraically  and  exponentially  for  Configs.  A  and  B 
respectively.  The  intermediate  grids  A2  and  B2  are  displayed  in  Fig.  4.2.  Some  of  the  salient 
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Table  4.1:  Flow  parameters 


Config. 

D 

Re  (x  106) 

Moo 

M'oo 

Q 

WM 

6 

x.hk 

Ythh 

*< 

3.0 

2.06 

5.25 

200.8 

20.6 

530 

-12.5 

-3.50 

-0.52 

B 

3.0 

1.53 

4.76 

222.6 

18.4 

530 

15.0 

Legend:  D  -  Diameter  of  cylinder  (in.) 

-  Upstream  Mach  Number 
-  Freestream  temperature  (R) 
Tw  -  Wall  temperature  (R) 


Re  -  Reynolds  number  (per  foot) 

-  Post-impinging  shock  Mach  Number 
Poo  -  Freestream  pressure  ( Ib/f 2) 

6  -  Flow  deflection  angle  (deg.) 


X,  Yahk  -  Ref-  location  of  shock  (in.) 


points  of  each  grid  are  presented  in  Table  4.2.  Guidelines  for  grid  resolution  are  taken  from  the 
work  of  Klopfer  and  Yee  [2]  who,  for  interfering  flows,  recommend  a  surface  cell  Reynolds  number 
( Rec)  of  roughly  3  for  heat-transfer  calculations  (fixed-wall  temperature)  and  of  the  order  of  10  for 
adiabatic  wall  conditions.  They  also  indicate  that  heat-transfer  rates  are  not  particularly  sensitive 
to  the  circumferential  spacing.  The  most  refined  grid  for  each  computation  satisfies  these  criteria. 
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Table  4.2:  Grid  details 


Legend:  IL  -  Points  in  £  direction  JL  -  Points  in  rj  direction 

Rec  -  Surface  mesh  Reynolds  number  A0  -  Angular  spacing  (deg) 
Subscripts:  av  -  average  min, max  -  minimum, maximum 
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Figure  4.1:  Geometry  and  boundary  conditions 


18 


Figure  4.2:  Grids  for  Configs.  A  and  B 
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5.  Boundary  Conditions  and  Numerical  Details 


With  reference  to  Fig.  4.1,  the  boundary  conditions  are  as  follows: 

•  ’’Inflow  boundary”  [BC]:  The  flow  vector  {/>,  pu,pv,pe}  is  specified  corresponding  to  the 
known  values,  either  upstream  or  downstream  of  the  impinging  shock. 

•  Solid  boundary  [AD]:  Since  this  is  a  solid  surface,  the  velocity  vector  and  the  normal  pressure 
gradient  are  assumed  zero  and  a  fixed  surface  temperature  is  specified  i.e., 

dv 

pv  =  0;  T  =  TW ;  ^  =  0  (5.1) 

where  the  subscript  w  refers  to  wall  conditions.  The  thermal  boundary  conditions  employed 
for  both  configurations  are  presented  in  Table  4.1. 

•  Outflow  boundaries  [AB,CD]:  The  flow  at  these  boundaries  is  assumed  predominantly  su¬ 
personic.  The  zero  gradient  extrapolation  condition  (d/d£  =  0)  is  applied. 

The  boundary  conditions  for  the  implicit  portion  of  the  algorithm  are  trivial  so  far  as  the 
freestream  is  concerned.  On  the  cylinder  surface  (j  =  1/2),  the  viscous  and  inviscid  terms  are 
treated  separately  to  account  for  the  split-flux  approach  by  appropriately  modifying  the  last  row  of 
the  block  tri-diagonal  system  represented  by  Eqn.  3.28  after  simplification  with  line  relaxation  [35]. 

In  all  computations  described,  convergence  is  determined  by  monitoring  several  quantities 
including  maximum  and  average  relative  change  in  the  solution  vector  over  a  flow  development 
time  of  roughly  three  characteristic  times  ( Tc  =  time  required  for  a  particle  to  traverse  the 
domain  at  upstream  conditions),  behavior  of  the  residual  and  surface  quantities  -  skin  friction 
(tw) ,  pressure  and  heat  transfer. 

With  default  compiler  vectorization,  the  code  executes  at  a  data  processing  rate  (DPR)  of 
roughly  1.17  x  10-3  secs/(iteration-grid  point)  on  a  CRAY-XMP  and  1.09  x  10“3  on  a  CRAY-2. 
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The  algorithm  is  very  stable  with  permissible  CFL  numbers  (based  on  inviscid  parameters)  above 
2000  and  5000  for  the  finer  and  coarser  grids  respectively.  This  permissible  CFL  is  approximately 
four  orders  of  magnitude  larger  than  that  achievable  with  the  explicit  MacCormack  unsplit  scheme. 
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6.  Results  and  Discussion 


6.1  Configuration  A 

The  computed  results  are  first  compared  with  experiment  and  the  computations  of  Thareja  et 
al.  [13].  Their  results  are  denoted  TSHMP  corresponding  to  the  first  initials  of  each  of  the 
authors.  They  computed  this  configuration  with  three  degrees  of  mesh  refinement  with  triangular 
elements.  In  the  following,  results  from  their  finest  grid,  consisting  of  5674  triangles  and  5792 
quadrilaterals,  are  compared  directly  with  those  from  the  finest  grid  utilized  in  this  research. 
The  current  algorithm  and  TSHMP  differ  significantly  in  approach  and  methodology.  A  precise 
comparison  of  numerical  efficiency  or  grid  sizes  employed  is  not  possible  since  unlike  the  current 
method,  Thareja  et  al.  utilize  unstructured  adapted  grids  in  their  approach. 

The  heat-transfer  comparison  is  presented  in  Fig.  6.1  in  which  are  plotted  the  normalized 
heat-transfer  values  against  angle  along  the  cylinder  (degrees)  measured  from  the  noninterfer¬ 
ing  stagnation  line  and  increasing  in  the  clockwise  direction.  For  normalization  purposes,  the 
approach  of  Thareja  et  al.  is  followed.  Although  the  experimental  noninterfering  stagnation 
heat-transfer  value  is  6l.7Btu/ ft 2  -  sec,  they  argue  that  since  the  theoretical  value  for  peak 
heat-transfer,  obtained  with  a  VSL  analysis  is  41.43,  the  difference  may  be  factored  out  by  uti¬ 
lizing  the  higher  (61.7)  value  to  normalize  experimental  results  and  the  lower  value  (41.43)  for 
computed  values.  The  results  (Fig.  6.1)  for  the  finest  grid  computed  indicate  good  agreement 
as  regards  both  peak  magnitude  and  peak  location  with  the  values  observed  by  TSHMP  with 
the  flux  corrected  transport  (FCT)  approach.  Peak  values  are  underpredicted  by  roughly  5%  by 
their  FCT  approach  and  by  roughly  6%  with  the  current  approach.  Their  basic  scheme  performs 
poorest.  Away  from  this  stagnation  value,  the  drop  in  heat-transfer  is  predicted  relatively  well  by 
both  methods  in  a  region  roughly  ±10°  around  the  stagnation  point.  Subsequently,  all  computed 
values  underpredict  heat-transfer  for  negative  9  values.  It  is  hypothesized  that  this  is  a  result 
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of  transition  to  turbulence.  A  small  bump  in  heat-transfer  values  is  observed  at  approximately 
-30  degrees  with  the  current  method.  Previous  research  for  the  Type  ///+  interaction  has  also 
displayed  similar  spatial  oscillations  downstream  of  the  stagnation  point  [12]. 

Fig.  6.2  compares  surface  pressure  values  normalized  with  the  noninterfering  stagnation  pres¬ 
sure.  Although  the  peak  pressure  location  precisely  matches  the  experimental  value  (to  within 
0.2  degrees),  the  value  itself  is  over  predicted  by  roughly  10%  on  the  finest  grid.  Further,  the  drop 
in  pressure  away  from  the  stagnation  point  is  more  rapid  with  the  current  method  than  observed 
in  experiment.  For  larger  negative  0  values,  significant  spatial  pressure  oscillation  is  observed 
qualitatively  similar  to  that  observed  for  the  heat-transfer  rates.  Although  TSHMP  results  do 
not  display  these  oscillations  at  least  two  other  methods  utilized  for  shock-shock  interactions  (van 
Leer  splitting  [12]  and  a  TVD  scheme  [2])  do  indeed  display  similar  behavior  for  Type  III  and 
Type  IV  flows.  For  each,  the  spatial  oscillations  occur  downstream  of  the  shear  layer  attachment 
described  below.  The  possibility  that  these  anomalies  are  slow-decaying  transients  is  currently 
under  investigation.  In  general  however,  the  comparison  of  the  current  method  with  experiment 
is  good  though  both  algorithms  directly  compared  tend  to  underpredict  pressures  away  from  the 
interaction. 

Fig.  6.3  displays  the  effect  of  grid  resolution  on  surface  pressure  and  heat-transfer  values.  Both 
peak  pressure  and  heat-transfer  values  increase  and  change  location  slightly  with  grid  refinement. 
For  negative  values  of  theta  and  away  from  the  stagnation  region,  the  better  resolved  grids  ( A2  and 
A3)  display  spatial  oscillations  as  observed  previously  for  causes  currently  unknown.  Furthermore, 
the  “bandwidth”  of  the  pressure  and  heat-transfer  profiles  also  reduces  as  several  features  of  the 
flow  (see  below)  are  captured  with  greater  clarity  in  the  finer  grids.  The  peak  values  and  locations 
for  both  pressure  and  heat-transfer  are  presented  in  Table  6.1.  The  computations  clearly  do  not 
display  the  desirable  grid  independence  and  further  investigation  is  warranted. 

The  ability  of  the  current  algorithm  to  extract  the  critical  flow  features  is  examined  by  com¬ 
parison  with  results  from  the  TSHMP  method.  For  the  purposes  of  brevity,  only  pressure  and 
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Mach  contours  are  discussed.  Fig.  6.4  displays  pressure  contours  for  the  three  computations  with 
the  current  method  and  the  TSHMP  method.  The  effect  of  grid  resolution  is,  as  anticipated,  to 
provide  clearer  definition  of  all  flow  features.  The  coarsest  grid  (Mesh  Al)  fails  to  resolve  the 
embedded  shock  (see  schematic  of  Fig.  1.2)  and  the  high-pressure  region  near  the  stagnation  point 
(because  of  the  terminating  jet  bow  shock)  where  the  supersonic  jet  impinges  on  the  surface  is 
diffuse.  The  largest  pressure  rise  in  the  system  occurs  across  the  distorted  bow  shock.  Mesh  A2 
clearly  captures  all  features  of  the  interaction  detectable  by  pressure  contours  although  it  is  evi¬ 
dent  from  the  surface  pressure  comparisons  presented  earlier  that  the  required  level  of  accuracy 
has  not  yet  been  reached.  The  highest  resolution  examined  in  this  research  (Mesh  3)  displays  the 
structure  of  the  embedded  shock  whose  orderly  deflection  by  a  terminating  compression  system  is 
clearly  exhibited.  This  compression  system  is  more  compact  with  the  current  method  than  with 
TSHMP.  An  examination  of  the  Mach  contours  (Fig.  6.5)  reiterates  the  conclusions  derived  from 
the  pressure  contours  in  so  far  as  the  shock  structure  is  concerned.  The  shear  layer  emanating 
from  the  intersection  point  of  the  impinging  shock  and  the  bow  shock  as  also  the  terminating 
jet  bow  shock  is  not  captured  at  all  with  the  coarsest  grid.  The  shear  layer  definition  emerges 
with  grid  refinement.  The  inviscid  shock  is  captured  within  at  most  two  grid  points  in  regions 
where  the  grid  is  aligned  with  the  shock  and  three  grid  points  elsewhere  as  anticipated  from  the 
choice  of  the  unmodified  first  order  Steger- Warming  splitting  in  the  vicinity  of  stronger  shock 

Table  6.1:  Peak  pressure  and  heat-transfer  values  -  Config.  A 


Case 

P/P0 

9  (deg.) 

Q/Qo 

9 

Al 

5.54 

-23.9 

5.5 

-23.9 

A2 

8.94 

-24.0 

-20.5 

A3 

9.70 

-24.5 

13.5 

-22.0 
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waves.  The  TSHMP  method  displays  crisper  bow  shock  resolution  that  may  be  at  least  partly 
attributable  to  grid  clustering  in  their  adaptive  method.  However,  the  finest  grid  in  the  current 
method  displays  better  resolution  of  the  shear  layer  structure  along  the  slip  and  body  surfaces. 
The  stand-off  distance  of  the  terminating  jet  bow  shock  is  slightly  higher  for  TSIIMP  than  with 
the  current  method. 

The  streamline  pattern  on  the  finest  grid  is  presented  in  Fig.  6.6.  A  stagnation  line  may  be 
clearly  defined  separating  particles  flowing  over  and  below  the  cylinder.  A  detailed  examination  of 
the  flow  reveals  the  existence  of  a  small  separated  region  of  recirculating  flow  above  the  horizontal 
axis  between  8  ~  21°  and  8  ~  34°  and  extending  over  roughly  30  grid  points  normal  to  the  wall. 
An  examination  of  the  surface  pressure  plot  (Fig.  6.2)  reveals  a  small  adverse  pressure  gradient  at 
8  ~  20°  and  most  likely  causes  separation  of  the  boundary  layer  developing  from  the  stagnation 
point.  Neither  of  the  coarser  grids  exhibit  this  separated  region. 

The  general  pattern  of  the  interaction  is  thus  clearly  exposed  in  the  current  computations  and 
may  be  classified  as  as  Type  IV.  The  impinging  shock  intersects  the  bow  shock  because  of  the 
cylinder  at  its  nearly  normal  portion.  A  supersonic  jet  is  formed  downstream  of  the  triple  point 
embedded  in  two  shear  layers.  The  size  of  the  supersonic  region  of  the  jet  may  be  deduced  from 
Fig.  6.7  which  shows  the  sonic  line  in  the  flow.  The  jet  then  terminates  in  a  jet  bow  shock  which 
in  turn  creates  a  region  of  stagnation  heating. 


6.2  Configuration  B 

The  results  obtained  with  the  current  scheme  for  a  Type  ///+  interaction  are  compared  with 
the  experimental  observations  of  Wieting  [39]  and  computations  with  the  van  Leer  flux-splitting 
scheme  with  and  without  turbulence  reported  by  Moon  and  Holt  [12].  Their  grid  is  closely 
comparable  to  the  intermediate  grid  utilized  in  this  research  (Grid  B2).  As  for  Config.  A,  surface 
quantities  are  first  examined  for  accuracy.  The  flow  is  then  analysed  for  critical  features. 
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Fig.  6.8  exhibits  surface  pressure  values  (normalized  by  the  noninterfering  stagnation  pressure) 
obtained  by  all  methods.  Experimental  values  remain  almost  constant  until  about  6  ~  25°  after 
which  they  display  a  steep  rise  upto  a  maximum  amplification  of  roughly  6.2  beyond  which 
they  drop  again  to  freestream  values.  The  laminar  computation  of  Moon  and  Holt  (abbreviated 
M  H  in  the  figures)  overpredicts  the  peak  value  by  about  100%  while  results  with  the  current 
method  utilizing  a  comparable  grid  (151  x  81  -  B2)  displays  overprediction  of  only  about  15%. 
Grid  refinement  with  the  current  method  overpredicts  pressure  with  a  maximum  amplification 
of  roughly  8.0.  The  best  peak  pressure  prediction  is  achieved  with  the  turbulent  computation  of 
Moon  and  Holt  where  the  peak  is  overpredicted  by  only  less  than  0.8  %.  All  methods  succesfully 
reproduce  the  extent  of  the  pressure  amplification  region  and  the  angular  position  of  the  pressure 
peak.  Both  computations  of  Moon  and  Holt  and  the  finest  grid  computed  with  the  current 
algorithm  exhibit  significant  spatial  waviness  in  the  circumferential  positions  between  -10°  and 
25°. 

Fig.  6.9  displays  a  comparison  of  heat-transfer  rate  prediction  where  the  ordinate  is  the  heat- 
transfer  normalized  by  experimentally  observed  stagnation  heat-transfer  values,  unlike  for  Con- 
fig.  A  where  results  from  a  VSL  analysis  were  utilized.  Away  from  the  peak  heating  region,  all 
laminar  computations  underpredict  heat-transfer  leading  to  the  speculation  that  the  boundary 
layer  on  either  side  of  the  stagnation  region  rapidly  transitions  into  turbulence.  The  computations 
of  Moon  and  Holt  (both  with  and  without  turbulence)  display  significant  spatial  oscillations  in 
the  region  immediately  adjacent  to  the  pressure  rise,  the  peak  experimental  heat-transfer  ampli¬ 
fication  value  is  approximately  6.6  occuring  at  37.9  degrees.  This  value  is  considerably  lower  than 
that  observed  previously  for  the  Type  IV  interaction.  Peak  heat-transfer  values  for  both  laminar 
and  turbulent  computations  of  Moon  and  Holt  underpredict  experiment  by  approximately  30%. 
On  comparable  grids,  (151  x  81  -  B2)  both  the  current  laminar  result  and  their  turbulent  com¬ 
putation  underpredict  the  peak  heat-transfer  although  the  current  method  performs  marginally 
better.  Witt  further  grid  refinement  (197  x  131  -  Case  B3),  peak  heat-transfer  values  overpredict 
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experimental  values  by  10%.  The  error  in  location  of  peak  heat-transfer  (roughly  3.0  degrees) 
may  be  attributed  to  small  errors  in  the  location  of  the  impinging  shock  which  has  previously 
been  observed  to  be  very  critical  [13].  Peak  values  and  their  locations  are  presented  in  Table  6.2. 

One  anomaly  evident  from  Fig.  6.9  is  the  presence  of  two  peaks,  of  equal  magnitude  for  case 
B2  and  unequal  magnitude  for  B3.  Since  the  turbulent  computation  of  Moon  and  Holt  also 
displays  this  behavior  to  some  extent,  this  phenomenon  was  examined  in  further  detail.  The  peak 
heat-transfer  rates  occur  at  the  point  where  the  jet  (Type  IV)  or  shear  layer  (Type  III)  impinges 
the  surface  of  the  model  and  depends  not  only  on  the  peak  pressure  generated  by  the  jet  but 
also  the  width  of  the  jet,  the  angle  at  which  the  jet  impinges  the  surface  and  whether  the  jet 
shear  layers  are  laminar  or  turbulent  [5].  A  superposition  of  Figs.  6.8  and  6.9  indicates  that  the 
single  pressure  peak  occurs  on  the  stagnation  line  which  lies  between  the  two  heat-transfer  peaks. 
The  relative  location  of  the  peaks  is  displayed  in  Fig.  6.10  which  shows  flow  streamlines  near  the 
stagnation  point  on  an  enlarged  scale.  One  possible  physical  explanation  hypothesizes  significant 
energy  transfer  between  the  flow  across  the  slipstream  emanating  from  the  triple  point.  As  a 
result,  the  enthalpy  of  the  flow  away  from  the  stagnation  line  is  concievably  augmented  which 
might  lead  to  higher  heat-transfer  at  points  A  and  B  in  Fig.  6.10.  Two  peaks  could  be  possible 
therefore  depending  upon  the  relative  location  of  the  shear  layer  and  the  stagnation  line.  Not 
much  information  can  be  derived  from  the  absence  of  a  pair  of  peaks  for  the  coarsest  grid  employed 

Table  6.2:  Peak  pressure  and  heat-transfer  values  -  Config.  B 


Case 

P/Po 

9  (deg.) 

Q/Qo 

9 

B1 

5.49 

46.0 

4.2 

43.6 

B2 

7.69 

45.7 

5.0 

43.9 

B3 

10.18 

45.0 

44.0 
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(B1  -  not  shown)  since  poor  resolution  smears  out  most  of  the  flow  structure. 

It  is  clear  that  the  above  mechanism  for  creation  of  the  second  peak  could  have  a  numerical 
origin.  To  investigate  this  possibility,  two  numerical  drawbacks  of  the  current  method  were 
examined.  First,  since  the  double  peaks  occur  in  the  vicinity  of  the  stagnation  point,  the  property 
of  flux  discontinuity  exhibited  by  the  current  algorithm  was  investigated.  The  eigenvalues  (see 
Eqn.  3.18)  were  smoothed  by  adding  terms  of  the  type  c  x  where  £  is  appropriately  either  f 
(for  F)  or  77  (for  G)  and  A  is  the  eigenvalue  of  the  flux  vector  under  consideration.  The  value 
of  c  was  fixed  at  0.01.  The  effect  of  this  smoothing  on  the  surface  properties  and  in  fact  the 
entire  solution  was  negligible  and  is  not  shown.  The  second  numerical  possibility  concerns  the  LR 
correction  described  in  Section  3.  As  explained  there,  the  effect  of  the  MacCormack  and  Candler 
modification  was  to  eliminate  a  first  order  damping  effect  at  the  expense  of  a  second  order  artificial 
pressure  gradient  whose  value  was  such  (Eqn.  3.24)  that  it  achieved  significant  values  only  under 
boundary  layer  assumptions.  The  LR  correction  was  proposed  as  a  solution,  and  in  fact,  all 
the  computations  described  with  the  current  method  thus  far  utilized  this  correction  in  a  linear 
fashion  but  only  for  the  surface- normal  77  direction.  A  shear  layer  such  as  that  emanating  from 
a  triple-point  and  impinging  on  a  surface  clearly  displays  significant  gradients  of  the  velocities 
in  the  circumferential  £  direction  and  correspondingly  a  fictitious  pressure  gradient  of  unknown 
magnitude  may  be  expected  to  exist  across  it.  A  small  region  extending  0.002  radii  from  the 
boundary  and  spanning  a  region  5°  on  either  side  of  the  stagnation  point  was  therefore  identified 
and  a  correction  similar  to  the  LR  correction  was  applied  also  in  the  f  direction.  The  results 
are  shown  in  Fig.  6.11  for  both  Cases  B2  and  B3.  The  anomaly  near  the  stagnation  point 
is  eliminated.  For  the  medium  grid  B2,  the  maximum  heat-transfer  increases  by  about  10% 
(underpredicts  maximum  experimental  values  by  16%)  while  for  the  finer  grid  B3,  the  effect  is 
much  more  pronounced  with  peak  heat  rates  increasing  nearly  25%  (overpredicting  maximum 
experimental  heat-transfer  by  45%).  The  location  of  the  peaks  is  2.5°  in  error.  The  surface 
pressure  is  not  significantly  influenced  and  is  therefore  not  shown. 


28 


Further  comparison  with  the  turbulent  computation  of  Moon  and  Holt  are  made  by  examining 
the  flow  structure.  The  equivalent  Grid  B2  is  utilized  for  this  purpose.  Fig.  6.12  displays  static 
pressure  contours  in  the  stagnation  region.  The  overall  comparison  is  good.  The  shock  structure 
is  identical  (to  within  one  grid  point)  although  the  smearing  of  the  shock  is  modestly  greater  with 
the  current  method  since  this  algorithm  reverts  to  the  original  Steger  Warming  formulation  at 
shock  waves.  The  reflected  shock  manifests  as  a  compression  fan  near  the  stagnation  point.  This 
fan  is  larger  for  the  current  method. 

The  Mach  contour  plots  are  displayed  in  Fig.  6.13  on  roughly  identical  scales.  The  flow  is 
generally  classified  as  a  type  ///+  since  it  possesses  features  representing  transition  from  type 
III  to  IV.  The  impinging  shock  interacts  with  the  bow  shock  causing  displacement  of  the  latter 
with  the  creation  of  a  transmitted  shock.  A  shear  layer  emanates  from  the  “triple”  point.  The 
flow  in  the  shear  layer  is  decelerated  as  it  approaches  the  surface  becoming  subsonic  with  a  near 
normal  shock  as  observed  in  Fig.  6.14  which  displays  the  sonic  line. 

Overall,  van  Leer’s  splitting  provides  a  sharper  representation  of  the  shocks  than  the  current 
scheme  which  is  essentially  a  Steger  Warming  scheme  near  shock  waves.  This  excessive  difFusi vity 
of  the  Steger  Warming  scheme  was  observed  previously  by  Anderson  et  al.  [32].  The  reflected 
shock  as  well  as  the  shear  layer  emanating  from  the  point  of  intersection  of  the  reflected  shock 
and  the  displaced  bow  shock  are  clearly  resolved  with  the  current  algorithm  although  over  more 
zones  than  with  van  Leer’s  splitting.  No  over  or  undershoots  in  Mach  number  are  observed  in  the 
current  results. 

Flow  streamlines  obtained  with  the  two  methods  are  displayed  in  Fig.  6.15.  Streams  deflected 
through  the  distorted  bow  shock  are  turned  toward  the  body  surface  but,  unlike  the  Type  IV 
interaction,  do  not  impinge  directly  upon  the  body  surface.  Moon  and  Holt  indicate  that  the 
features  separating  a  Type  III+  flow  from  a  Type  III  flow  are:  (i)  infant  stage  of  the  supersonic 
jet,  (ii)  terminating  bow  shock,  and  (iii)  boundary  layer  separation.  The  infant  stage  of  the 
supersonic  (feature  (i))  jet  is  observed  in  the  flow  above  the  stagnation  line  which  first  decelerates 
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through  a  strong  shock  (feature  (ii)).  For  this  interaction,  the  shock  is  nearly  perpendicular  to  the 
body  surface  because  of  the  particular  angle  between  the  shear  layer  flow  and  the  body  surface. 
Subsequently,  the  flow  rapidly  accelerates  to  supersonic  and  may  be  considered  to  be  the  inception 
of  a  supersonic  jet  that  would  exist  if  the  interaction  were  Type  IV.  This  may  happen  if,  for 
example,  the  shear  layer  angle  changed  modestly  in  the  direction  of  the  body  so  as  to  graze  the 
cylinder  surface.  Regarding  feature  (iii),  no  boundary  layer  separation  was  observed,  for  any  of 
the  computations  performed  for  this  configuration. 


30 


Figure  6.1:  Surface  heat-transfer  comparison  -  Config.  A 
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Figure  6.2:  Surface  pressure  comparison  -  Config.  A 
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(a)  Surface  Pressure 


(b)  Surface  heat  transfer 


Figure  6.3:  Effect  of  grid  resolution  -  Config.  A 
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Figure  6.4:  Comparison  of  pressure  contours  -  Config.  A  (..over) 
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Figure  6.5:  Comparison  of  Mach  contours  -  Config.  A  (..over) 
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Figure  6.5:  Comparison  of  Mach  contours  -  Config.  A  (concluded) 
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Figure  6.6:  Streamlines  -  Config.  A  (Grid  A3) 
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Figure  6.8:  Comparison  of  surface  pressure  -  Config.  B 


Figure  6.9:  Comparison  of  surface  heat-transfer  -  Config.  B 
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Figure  6.12:  Comparison  of  pressure  contours  -  Config.  B 
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Figure  6.13:  Comparison  of  Mach  contours  -  Config.  B 
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Figure  6.15:  Streamline  comparison  -  Config.  B 
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7.  Conclusions 


Viscous  shock-on-shock  interactions  are  examined  with  a  modified  Steger  Warming  scheme  dis¬ 
playing  low  numerical  dissipation  in  boundary  layers.  For  a  Type  IV  interaction,  peak  heat 
transfer  and  pressure  results  are  in  relatively  good  agreement  with  experiment.  Modest  discrep¬ 
ancies  away  from  the  stagnation  point,  we  believe,  are  due  to  the  existence  of  boundary  layer 
transition  ^  turbulence  which  was  not  modeled  in  this  research.  The  computed  shock  and  shear 
layer  patterns  compare  well  with  those  obtained  with  other  methods.  A  small  region  of  separated 
flow  is  observed  in  the  finest  grid  computed.  For  a  Type  III+  interaction,  the  algorithm  displays 
spatial  heat-transfer  oscillations  caused  most  likely  by  fictitious  pressure  gradients  in  the  region 
of  peak  heating.  These  gradients  are  successfully  eliminated  by  appropriate  local  corrections. 
Overall,  grid  refinement  leads  to  overprediction  of  peak  pressures  and  heat-transfer  rates.  Thus, 
on  the  meshes  examined,  grid  independence  is  not  observed  and  further  investigation  is  necessary. 
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Nomenclature 

Jacobian  of  flux  vector  G  with  respect  to  U 
generic  constant 

specific  heat  at  constant  pressure 
specific  heat  at  constant  volume 
diameter 

Data  Processing  Rate 
total  energy 
internal  energy 
flux  vectors 

Flux  Corrected  Transport  Method 

Points  in  £  and  rj  directions  respectively 

Jacobian  of  transformation 

Last  Row  correction 

Mach  number 

iteration  number 

pressure 

heat  transfer  component;  similarity  transformation  matrix 
Gas  constant  for  air 
Reynolds  number 
time 

Temperature;  time 

Thareja,  Stewart,  Hassan,  Morgan  and  Peraire 
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TVD  Total  Variation  Diminishing  Schemes 
u  Cartesian  velocity  in  x  direction 

U  solution  vector 

U  contravariant  velocity  vector 

v  Cartesian  velocity  in  y  direction 

v  Cartesian  velocity  vector 

V  contravariant  velocity  vector 

VSL  Viscous  Shear  Layer  Analysis 

x,y  physical  coordinates 

a  kinetic  energy 

6  change;  angle 

A  switch  threshold;  change 

d  partial  derivative 

rj  transformed  coordinate 

7  ratio  of  specific  heats 

A  second  coefficient  of  viscosity 

p.  molecular  viscosity 

<t>  constant  in  limiter  function 

p  density 

t  Cartesian  stress  tensor 

9  angle 
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£  transformed  variable 


Subscripts 
av  average 
c  characteristic 
i  component  i;  internal 
j  component  j;  face  j 
max  maximum  value 
min  minimum  value 
w  wall 

x,  y  Cartesian  coordinate  directions 
77  derivative  with  respect  to  r/ 

{  derivative  with  respect  to  f 

00  freestream  value 
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